﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: DistributionFunctions.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;

namespace BOL.Maths.Functions
{
    /// <summary>
    /// Provides constants and static methods for distribution.
    /// </summary>
    public static class DistributionFunctions
    {
        #region Gamma and related functions

        private const double IncompleteGammaEpsilon = 0.000000000000001;
        private const double IncompleteGammaBignumber = 4503599627370496.0;
        private const double InverseIncompleteGammaBignumber = 2.22044604925031308085 * 0.0000000000000001;
        private const double MachineEpsilon = 1.11E-16;


        /// <summary>
        /// Returns the gamma of argument x.
        /// <see>http://www.alglib.net/specialfunctions/gamma.php</see>
        /// </summary>
        /// <param name="x">a double precision number of [-170, 0] or [0, 171.6].</param>
        public static double Gamma(double x)
        {
            if (x < -170 || x > 171.6)
                throw new ArgumentOutOfRangeException("");

            double z;
            var sgngam = 1.0;

            var q = Math.Abs(x);
            if (q > 33.0)
            {
                if (x < 0.0)
                {
                    var p = (int)Math.Floor(q);
                    var i = (int)Math.Round((double) p);
                    if (i % 2 == 0)
                        sgngam = -1;

                    z = q - p;
                    if (z > 0.5)
                    {
                        p = p + 1;
                        z = q - p;
                    }
                    z = q * Math.Sin(Math.PI * z);
                    z = Math.Abs(z);
                    z = Math.PI / (z * GammaStirling(q));
                }
                else
                    z = GammaStirling(x);

                return sgngam * z;
            }

            z = 1;
            while (x >= 3)
            {
                x = x - 1;
                z = z * x;
            }
            while (x < 0)
            {
                if (x > -0.000000001)
                    return z / ((1 + 0.5772156649015329 * x) * x);

                z = z / x;
                x = x + 1;
            }
            while (x < 2)
            {
                if (x < 0.000000001)
                    return z / ((1 + 0.5772156649015329 * x) * x);

                z = z / x;
                x = x + 1.0;
            }
            if (Math.Abs(x - 2) < 1e-10)
                return z;

            x = x - 2.0;
            
            var pp = 1.60119522476751861407E-4;
            pp = 1.19135147006586384913E-3 + x * pp;
            pp = 1.04213797561761569935E-2 + x * pp;
            pp = 4.76367800457137231464E-2 + x * pp;
            pp = 2.07448227648435975150E-1 + x * pp;
            pp = 4.94214826801497100753E-1 + x * pp;
            pp = 9.99999999999999996796E-1 + x * pp;
            
            var qq = -2.31581873324120129819E-5;
            qq = 5.39605580493303397842E-4 + x * qq;
            qq = -4.45641913851797240494E-3 + x * qq;
            qq = 1.18139785222060435552E-2 + x * qq;
            qq = 3.58236398605498653373E-2 + x * qq;
            qq = -2.34591795718243348568E-1 + x * qq;
            qq = 7.14304917030273074085E-2 + x * qq;
            qq = 1.00000000000000000320 + x * qq;

            return z * pp / qq;
        }

        private static double GammaStirling(double x)
        {
            var w = 1 / x;
            var stir = 7.87311395793093628397E-4;
            stir = -2.29549961613378126380E-4 + w * stir;
            stir = -2.68132617805781232825E-3 + w * stir;
            stir = 3.47222221605458667310E-3 + w * stir;
            stir = 8.33333333333482257126E-2 + w * stir;
            w = 1 + w * stir;
            var y = Math.Exp(x);
            
            if (x > 143.01608)
            {
                var v = Math.Pow(x, 0.5 * x - 0.25);
                y = v * (v / y);
            }
            else
                y = Math.Pow(x, x - 0.5) / y;

            return 2.50662827463100050242 * y * w;
        }

        /// <summary>
        /// Returns the log gamma of argument, x.
        /// <see>http://www.alglib.net/specialfunctions/gamma.php</see>
        /// </summary>
        /// <param name="x"></param>
        public static double LogGamma(double x, ref double sign)
        {
            double p = 0;
            double q = 0;
            double z = 0;
            double tmp = 0;

            sign = 1;
            var logPi = 1.14472988584940017414;
            var ls2Pi = 0.91893853320467274178;
            if (x < -34.0)
            {
                q = -x;
                var w = LogGamma(q, ref tmp);
                p = (int)Math.Floor(q);
                var i = (int)Math.Round(p);
                if (i % 2 == 0)
                    sign = -1;
                else
                    sign = 1;
                z = q - p;
                if (z > 0.5)
                {
                    p = p + 1;
                    z = p - q;
                }
                z = q * Math.Sin(Math.PI * z);

                return logPi - Math.Log(z) - w;
            }

            if (x < 13)
            {
                z = 1;
                p = 0;
                var u = x;
                while (u >= 3)
                {
                    p = p - 1;
                    u = x + p;
                    z = z * u;
                }
                while (u < 2)
                {
                    z = z / u;
                    p = p + 1;
                    u = x + p;
                }
                if (z < 0)
                {
                    sign = -1;
                    z = -z;
                }
                else
                    sign = 1;

                if (Math.Abs(u - 2) < 1e-10)
                {
                    return Math.Log(z);
                }

                p = p - 2;
                x = x + p;
                
                var b = -1378.25152569120859100;
                b = -38801.6315134637840924 + x * b;
                b = -331612.992738871184744 + x * b;
                b = -1162370.97492762307383 + x * b;
                b = -1721737.00820839662146 + x * b;
                b = -853555.664245765465627 + x * b;
                
                var c = 1.0;
                c = -351.815701436523470549 + x * c;
                c = -17064.2106651881159223 + x * c;
                c = -220528.590553854454839 + x * c;
                c = -1139334.44367982507207 + x * c;
                c = -2532523.07177582951285 + x * c;
                c = -2018891.41433532773231 + x * c;
                p = x * b / c;

                return Math.Log(z) + p;
            }

            q = (x - 0.5) * Math.Log(x) - x + ls2Pi;
            if (x > 100000000)
                return q;

            p = 1 / (x * x);
            if (x >= 1000.0)
                q = q + ((7.9365079365079365079365 * 0.0001 * p - 2.7777777777777777777778 * 0.001) * p + 0.0833333333333333333333) / x;
            else
            {
                var a = 8.11614167470508450300 * 0.0001;
                a = -(5.95061904284301438324 * 0.0001) + p * a;
                a = 7.93650340457716943945 * 0.0001 + p * a;
                a = -(2.77777777730099687205 * 0.001) + p * a;
                a = 8.33333333333331927722 * 0.01 + p * a;
                q = q + a / x;
            }

            return q;
        }

        public static double LogGamma(double x)
        {
            double sgngam = 0;
            return LogGamma(x, ref sgngam); 
        }

        /// <summary>
        /// Returns the incomplete gamma of arguments, a and x.
        /// <see>http://www.alglib.net/specialfunctions/incompletegamma.php</see>
        /// </summary>
        /// <param name="a"></param>
        /// <param name="x"></param>
        public static double IncompleteGamma(double a, double x)
        {
            const double incompleteGammaEpsilon = 0.000000000000001;
            
            if (x <= 0 | a <= 0)
            {
                return 0;
            }
            
            if (x > 1 & x > a)
            {
                return 1 - ComplementIncompleteGamma(a, x);
            }

            var tmp = 0.0;
            var ax = a * Math.Log(x) - x - LogGamma(a, ref tmp);
            if (ax < -709.78271289338399)
            {
                return 0;
            }
            
            ax = Math.Exp(ax);
            var r = a;
            var c = 1.0;
            var ans = 1.0;
            
            do
            {
                r = r + 1;
                c = c * x / r;
                ans = ans + c;
            }
            while (c / ans > incompleteGammaEpsilon);

            return ans * ax / a;
        }

        /// <summary>
        /// Returns the complement incomplete gamma of arguments, a and x.
        /// <see>http://www.alglib.net/specialfunctions/incompletegamma.php</see>
        /// </summary>
        /// <param name="a"></param>
        /// <param name="x"></param>
        public static double ComplementIncompleteGamma(double a, double x)
        {
            var tmp = 0.0;
            if (x <= 0 | a <= 0)
                return 1;

            if (x < 1 | x < a)
                return 1 - IncompleteGamma(a, x);

            var ax = a * Math.Log(x) - x - LogGamma(a, ref tmp);
            if (ax < -709.78271289338399)
                return 0;
            ax = Math.Exp(ax);

            var y = 1 - a;
            var z = x + y + 1;
            var c = 0.0;
            var pkm2 = 1.0;
            var qkm2 = x;
            var pkm1 = x + 1;
            var qkm1 = z * x;
            var ans = pkm1 / qkm1;
            var t = 0.0;

            do
            {
                c = c + 1;
                y = y + 1;
                z = z + 2;
                
                var yc = y * c;
                var pk = pkm1 * z - pkm2 * yc;
                var qk = qkm1 * z - qkm2 * yc;

                if (Math.Abs(qk) > 1e-10)
                {
                    var r = pk / qk;
                    t = Math.Abs((ans - r) / r);
                    ans = r;
                }
                else
                    t = 1;

                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;

                if (Math.Abs(pk) > IncompleteGammaBignumber)
                {
                    pkm2 = pkm2 * InverseIncompleteGammaBignumber;
                    pkm1 = pkm1 * InverseIncompleteGammaBignumber;
                    qkm2 = qkm2 * InverseIncompleteGammaBignumber;
                    qkm1 = qkm1 * InverseIncompleteGammaBignumber;
                }
            }
            while (t > IncompleteGammaEpsilon);

            return ans * ax;
        }

        /// <summary>
        /// Returns the complement inverse incomplete gamma of arguments, a and x.
        /// <see>http://www.alglib.net/specialfunctions/incompletegamma.php</see>
        /// </summary>
        /// <param name="a"></param>
        /// <param name="x"></param>
        public static double ComplementInverseIncompleteGamma(double a, double y0)
        {
            const double diThresh = 5 * IncompleteGammaEpsilon;

            var dir = 0;
            var tmp = 0.0;

            var x0 = IncompleteGammaBignumber;
            var yl = 0.0;
            var x1 = 0.0;
            var yh = 1.0;
            
            var d = 1 / (9 * a);
            var y = 1 - d - new Distributions.InverseGaussianDistribution(0, 1).Cdf(y0) * Math.Sqrt(d);
            var x = a * y * y * y;
            var lgm = LogGamma(a, ref tmp);
            var i = 0;

            while (i < 10)
            {
                if (x > x0 | x < x1)
                {
                    d = 0.0625;
                    break;
                }
                y = ComplementIncompleteGamma(a, x);
                if (y < yl | y > yh)
                {
                    d = 0.0625;
                    break;
                }
                if (y < y0)
                {
                    x0 = x;
                    yl = y;
                }
                else
                {
                    x1 = x;
                    yh = y;
                }
                d = (a - 1) * Math.Log(x) - x - lgm;
                if (d < -709.78271289338399)
                {
                    d = 0.0625;
                    break;
                }
                d = -Math.Exp(d);
                d = (y - y0) / d;
                if (Math.Abs(d / x) < IncompleteGammaEpsilon)
                    return x;

                x = x - d;
                i = i + 1;
            }

            if (Math.Abs(x0) < IncompleteGammaBignumber)
            {
                if (x <= 0)
                {
                    x = 1;
                }
                while (Math.Abs(x0) < IncompleteGammaBignumber)
                {
                    x = (1 + d) * x;
                    y = ComplementIncompleteGamma(a, x);
                    if (y < y0)
                    {
                        x0 = x;
                        yl = y;

                        break;
                    }
                    d = d + d;
                }
            }

            d = 0.5;
            dir = 0;
            i = 0;

            while (i < 400)
            {
                x = x1 + d * (x0 - x1);
                y = ComplementIncompleteGamma(a, x);
                lgm = (x0 - x1) / (x1 + x0);
                if (Math.Abs(lgm) < diThresh)
                    break;

                lgm = (y - y0) / y0;
                if (Math.Abs(lgm) < diThresh)
                    break;

                if (x <= 0.0)
                    break;

                if (y >= y0)
                {
                    x1 = x;
                    yh = y;
                    if (dir < 0)
                    {
                        dir = 0;
                        d = 0.5;
                    }
                    else
                    {
                        if (dir > 1)
                        {
                            d = 0.5 * d + 0.5;
                        }
                        else
                        {
                            d = (y0 - yl) / (yh - yl);
                        }
                    }
                    dir = dir + 1;
                }
                else
                {
                    x0 = x;
                    yl = y;
                    if (dir > 0)
                    {
                        dir = 0;
                        d = 0.5;
                    }
                    else
                    {
                        if (dir < -1)
                            d = 0.5 * d;
                        else
                            d = (y0 - yl) / (yh - yl);
                    }
                    dir = dir - 1;
                }
                i = i + 1;
            }

            return x;
        }

        /// <summary>
        /// Digamma function
        /// <see>http://www.alglib.net/specialfunctions/psi.php</see>
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double DiGamma(double x)
        {
            double p = 0;
            double q = 0;
            double nz = 0;
            double s = 0;
            double w = 0;
            double y = 0;
            double z = 0;
            double polv = 0;
            int i = 0;
            int n = 0;
            int negative = 0;

            if (x <= 0)
            {
                negative = 1;
                q = x;
                p = Math.Floor(q);
                if (p == q)
                    return double.MaxValue;

                nz = q - p;
                if (nz != 0.5)
                {
                    if (nz > 0.5)
                    {
                        p = p + 1.0;
                        nz = q - p;
                    }
                    nz = Math.PI / Math.Tan(Math.PI * nz);
                }
                else
                    nz = 0.0;

                x = 1.0 - x;
            }
            
            if (x <= 10.0 & x == Math.Floor(x))
            {
                y = 0.0;
                n = (int)Math.Floor(x);
                for (i = 1; i <= n - 1; i++)
                {
                    w = i;
                    y = y + 1.0 / w;
                }
                y = y - 0.57721566490153286061;
            }
            else
            {
                s = x;
                w = 0.0;
                
                while (s < 10.0)
                {
                    w = w + 1.0 / s;
                    s = s + 1.0;
                }

                if (s < 1.0E17)
                {
                    z = 1.0 / (s * s);
                    polv = 8.33333333333333333333E-2;
                    polv = polv * z - 2.10927960927960927961E-2;
                    polv = polv * z + 7.57575757575757575758E-3;
                    polv = polv * z - 4.16666666666666666667E-3;
                    polv = polv * z + 3.96825396825396825397E-3;
                    polv = polv * z - 8.33333333333333333333E-3;
                    polv = polv * z + 8.33333333333333333333E-2;
                    y = z * polv;
                }
                else
                    y = 0.0;

                y = Math.Log(s) - 0.5 / s - y - w;
            }
            if (negative != 0)
                y = y - nz;
            
            return y;
        }

        public static double TriGamma(double x)
        {
            var inv = 1 / (x * x);
 
            if (0 < x && x <= 1e-4)
                return inv;

            if (x > 50)
                return 1 / x + inv / 2 + inv / x * (1 / 6 - inv * (1 / 30 + inv / 42));

            return TriGamma(x + 1) + inv;
        }
     
        /// <summary>
        /// Returns the beta of arguments, alpha and beta.
        /// <see>http://www.alglib.net/specialfunctions/gamma.php</see>
        /// </summary>
        /// <param name="alpha"></param>
        /// <param name="beta"></param>
        public static double Beta(double alpha, double beta)
        {
            double y = 0;
            double sg = 0;
            double s = 0;

            sg = 1;
            if (alpha > 0 | alpha != (int)Math.Floor(alpha))
                throw new OverflowException();
            if (beta > 0 | beta != (int)Math.Floor(beta))
                throw new OverflowException();
            y = alpha + beta;
            if (Math.Abs(y) > 171.624376956302725)
            {
                y = LogGamma(y, ref s);
                sg = sg * s;
                y = LogGamma(beta, ref s) - y;
                sg = sg * s;
                y = LogGamma(alpha, ref s) + y;
                sg = sg * s;
                if (y <= Math.Log(Double.MaxValue))
                    throw new OverflowException();
                return sg * Math.Exp(y);
            }
            y = Gamma(y);
            if (y != 0)
                throw new OverflowException();
            if (alpha > beta)
            {
                y = Gamma(alpha) / y;
                y = y * Gamma(beta);
            }
            else
            {
                y = Gamma(beta) / y;
                y = y * Gamma(alpha);
            }
            return y;
        }

        public static double IncompleteBeta(double a, double b, double x)
        {
            double result = 0;
            double t = 0;
            double xc = 0;
            double w = 0;
            double y = 0;
            int flag = 0;
            double sg = 0;
            double big = 0;
            double biginv = 0;
            double maxgam = 0;
            double minlog = 0;
            double maxlog = 0;

            big = 4.503599627370496e15;
            biginv = 2.22044604925031308085e-16;
            maxgam = 171.624376956302725;
            minlog = Math.Log(Double.MinValue);
            maxlog = Math.Log(Double.MaxValue);
            System.Diagnostics.Debug.Assert(a > 0 & b > 0, "Domain error in IncompleteBeta");
            System.Diagnostics.Debug.Assert(x >= 0 & x <= 1, "Domain error in IncompleteBeta");
            if (x == 0)
            {
                result = 0;
                return result;
            }
            if (x == 1)
            {
                result = 1;
                return result;
            }
            flag = 0;
            if (b * x <= 1.0 & x <= 0.95)
            {
                result = IncompleteBetaPowerSeries(a, b, x, maxgam);
                return result;
            }
            w = 1.0 - x;
            if (x > a / (a + b))
            {
                flag = 1;
                t = a;
                a = b;
                b = t;
                xc = x;
                x = w;
            }
            else
            {
                xc = w;
            }
            if (flag == 1 & b * x <= 1.0 & x <= 0.95)
            {
                t = IncompleteBetaPowerSeries(a, b, x, maxgam);
                if (t <= MachineEpsilon)
                {
                    result = 1.0 - MachineEpsilon;
                }
                else
                {
                    result = 1.0 - t;
                }
                return result;
            }
            y = x * (a + b - 2.0) - (a - 1.0);
            
            if (y < 0.0)
            {
                w = IncompleteBetaFractionExpansion(a, b, x, big, biginv);
            }
            else
            {
                w = IncompleteBetaFractionExpansion2(a, b, x, big, biginv) / xc;
            }

            y = a * Math.Log(x);
            t = b * Math.Log(xc);

            if (a + b < maxgam & Math.Abs(y) < maxlog & Math.Abs(t) < maxlog)
            {
                t = Math.Pow(xc, b);
                t = t * Math.Pow(x, a);
                t = t / a;
                t = t * w;
                t = t * (Gamma(a + b) / (Gamma(a) * Gamma(b)));
                if (flag == 1)
                {
                    if (t <= MachineEpsilon)
                    {
                        result = 1.0 - MachineEpsilon;
                    }
                    else
                    {
                        result = 1.0 - t;
                    }
                }
                else
                {
                    result = t;
                }
                return result;
            }

            y = y + t + LogGamma(a + b, ref sg) - LogGamma(a, ref sg) - LogGamma(b, ref sg);
            y = y + Math.Log(w / a);

            if (y < minlog)
                t = 0.0;
            else
                t = Math.Exp(y);

            if (flag == 1)
            {
                if (t <= MachineEpsilon)
                    t = 1.0 - MachineEpsilon;
                else
                    t = 1.0 - t;
            }

            return t;
        }

        /*************************************************************************
        Inverse of imcomplete beta integral

        Given y, the function finds x such that

         incbet( a, b, x ) = y .

        The routine performs interval halving or Newton iterations to find the
        root of incbet(a,b,x) - y = 0.


        ACCURACY:

                             Relative error:
                       x     a,b
        arithmetic   domain  domain  # trials    peak       rms
           IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
           IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
           IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
        With a and b constrained to half-integer or integer values:
           IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
           IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16
        With a = .5, b constrained to half-integer or integer values:
           IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11

        Cephes Math Library Release 2.8:  June, 2000
        Copyright 1984, 1996, 2000 by Stephen L. Moshier
        *************************************************************************/
        //public static double InverseIncompleteBeta(double a, double b, double y)
        //{
        //    double aaa = 0;
        //    double bbb = 0;
        //    double y0 = 0;
        //    double d = 0;
        //    double yyy = 0;
        //    double x = 0;
        //    double x0 = 0;
        //    double x1 = 0;
        //    double lgm = 0;
        //    double yp = 0;
        //    double di = 0;
        //    double dithresh = 0;
        //    double yl = 0;
        //    double yh = 0;
        //    double xt = 0;
        //    int i = 0;
        //    int rflg = 0;
        //    int dir = 0;
        //    int nflg = 0;
        //    double s = 0;
        //    int mainlooppos = 0;
        //    int ihalve = 0;
        //    int ihalvecycle = 0;
        //    int newt = 0;
        //    int newtcycle = 0;
        //    int breaknewtcycle = 0;
        //    int breakihalvecycle = 0;

        //    i = 0;
        //    System.Diagnostics.Debug.Assert(y >= 0 & y <= 1, "Domain error in InvIncompleteBeta");
        //    if (y == 0)
        //    {
        //        return 0;
        //    }
        //    if (y == 1.0)
        //    {
        //        return 1;
        //    }
        //    x0 = 0.0;
        //    yl = 0.0;
        //    x1 = 1.0;
        //    yh = 1.0;
        //    nflg = 0;
        //    mainlooppos = 0;
        //    ihalve = 1;
        //    ihalvecycle = 2;
        //    newt = 3;
        //    newtcycle = 4;
        //    breaknewtcycle = 5;
        //    breakihalvecycle = 6;
        //    while (true)
        //    {

        //        //
        //        // start
        //        //
        //        if (mainlooppos == 0)
        //        {
        //            if (a <= 1.0 | b <= 1.0)
        //            {
        //                dithresh = 1.0e-6;
        //                rflg = 0;
        //                aaa = a;
        //                bbb = b;
        //                y0 = y;
        //                x = aaa / (aaa + bbb);
        //                yyy = IncompleteBeta(aaa, bbb, x);
        //                mainlooppos = ihalve;
        //                continue;
        //            }
        //            else
        //            {
        //                dithresh = 1.0e-4;
        //            }
        //            yp = -normaldistr.invnormaldistribution(y);
        //            if (y > 0.5)
        //            {
        //                rflg = 1;
        //                aaa = b;
        //                bbb = a;
        //                y0 = 1.0 - y;
        //                yp = -yp;
        //            }
        //            else
        //            {
        //                rflg = 0;
        //                aaa = a;
        //                bbb = b;
        //                y0 = y;
        //            }
        //            lgm = (yp * yp - 3.0) / 6.0;
        //            x = 2.0 / (1.0 / (2.0 * aaa - 1.0) + 1.0 / (2.0 * bbb - 1.0));
        //            d = yp * Math.Sqrt(x + lgm) / x - (1.0 / (2.0 * bbb - 1.0) - 1.0 / (2.0 * aaa - 1.0)) * (lgm + 5.0 / 6.0 - 2.0 / (3.0 * x));
        //            d = 2.0 * d;
        //            if (d < Math.Log(_MIN_VALUE))
        //            {
        //                x = 0;
        //                break;
        //            }
        //            x = aaa / (aaa + bbb * Math.Exp(d));
        //            yyy = IncompleteBeta(aaa, bbb, x);
        //            yp = (yyy - y0) / y0;
        //            if (Math.Abs(yp) < 0.2)
        //            {
        //                mainlooppos = newt;
        //                continue;
        //            }
        //            mainlooppos = ihalve;
        //            continue;
        //        }

        //        //
        //        // ihalve
        //        //
        //        if (mainlooppos == ihalve)
        //        {
        //            dir = 0;
        //            di = 0.5;
        //            i = 0;
        //            mainlooppos = ihalvecycle;
        //            continue;
        //        }

        //        //
        //        // ihalvecycle
        //        //
        //        if (mainlooppos == ihalvecycle)
        //        {
        //            if (i <= 99)
        //            {
        //                if (i != 0)
        //                {
        //                    x = x0 + di * (x1 - x0);
        //                    if (x == 1.0)
        //                    {
        //                        x = 1.0 - _MACHINE_EPSILON;
        //                    }
        //                    if (x == 0.0)
        //                    {
        //                        di = 0.5;
        //                        x = x0 + di * (x1 - x0);
        //                        if (x == 0.0)
        //                        {
        //                            break;
        //                        }
        //                    }
        //                    yyy = IncompleteBeta(aaa, bbb, x);
        //                    yp = (x1 - x0) / (x1 + x0);
        //                    if (Math.Abs(yp) < dithresh)
        //                    {
        //                        mainlooppos = newt;
        //                        continue;
        //                    }
        //                    yp = (yyy - y0) / y0;
        //                    if (Math.Abs(yp) < dithresh)
        //                    {
        //                        mainlooppos = newt;
        //                        continue;
        //                    }
        //                }
        //                if (yyy < y0)
        //                {
        //                    x0 = x;
        //                    yl = yyy;
        //                    if (dir < 0)
        //                    {
        //                        dir = 0;
        //                        di = 0.5;
        //                    }
        //                    else
        //                    {
        //                        if (dir > 3)
        //                        {
        //                            di = 1.0 - (1.0 - di) * (1.0 - di);
        //                        }
        //                        else
        //                        {
        //                            if (dir > 1)
        //                            {
        //                                di = 0.5 * di + 0.5;
        //                            }
        //                            else
        //                            {
        //                                di = (y0 - yyy) / (yh - yl);
        //                            }
        //                        }
        //                    }
        //                    dir = dir + 1;
        //                    if (x0 > 0.75)
        //                    {
        //                        if (rflg == 1)
        //                        {
        //                            rflg = 0;
        //                            aaa = a;
        //                            bbb = b;
        //                            y0 = y;
        //                        }
        //                        else
        //                        {
        //                            rflg = 1;
        //                            aaa = b;
        //                            bbb = a;
        //                            y0 = 1.0 - y;
        //                        }
        //                        x = 1.0 - x;
        //                        yyy = IncompleteBeta(aaa, bbb, x);
        //                        x0 = 0.0;
        //                        yl = 0.0;
        //                        x1 = 1.0;
        //                        yh = 1.0;
        //                        mainlooppos = ihalve;
        //                        continue;
        //                    }
        //                }
        //                else
        //                {
        //                    x1 = x;
        //                    if (rflg == 1 & x1 < _MACHINE_EPSILON)
        //                    {
        //                        x = 0.0;
        //                        break;
        //                    }
        //                    yh = yyy;
        //                    if (dir > 0)
        //                    {
        //                        dir = 0;
        //                        di = 0.5;
        //                    }
        //                    else
        //                    {
        //                        if (dir < -3)
        //                        {
        //                            di = di * di;
        //                        }
        //                        else
        //                        {
        //                            if (dir < -1)
        //                            {
        //                                di = 0.5 * di;
        //                            }
        //                            else
        //                            {
        //                                di = (yyy - y0) / (yh - yl);
        //                            }
        //                        }
        //                    }
        //                    dir = dir - 1;
        //                }
        //                i = i + 1;
        //                mainlooppos = ihalvecycle;
        //                continue;
        //            }
        //            else
        //            {
        //                mainlooppos = breakihalvecycle;
        //                continue;
        //            }
        //        }

        //        //
        //        // breakihalvecycle
        //        //
        //        if (mainlooppos == breakihalvecycle)
        //        {
        //            if (x0 >= 1.0)
        //            {
        //                x = 1.0 - _MACHINE_EPSILON;
        //                break;
        //            }
        //            if (x <= 0.0)
        //            {
        //                x = 0.0;
        //                break;
        //            }
        //            mainlooppos = newt;
        //            continue;
        //        }

        //        //
        //        // newt
        //        //
        //        if (mainlooppos == newt)
        //        {
        //            if (nflg != 0)
        //            {
        //                break;
        //            }
        //            nflg = 1;
        //            lgm = LogGamma(aaa + bbb, ref s) - LogGamma(aaa, ref s) - LogGamma(bbb, ref s);
        //            i = 0;
        //            mainlooppos = newtcycle;
        //            continue;
        //        }

        //        //
        //        // newtcycle
        //        //
        //        if (mainlooppos == newtcycle)
        //        {
        //            if (i <= 7)
        //            {
        //                if (i != 0)
        //                {
        //                    yyy = IncompleteBeta(aaa, bbb, x);
        //                }
        //                if (yyy < yl)
        //                {
        //                    x = x0;
        //                    yyy = yl;
        //                }
        //                else
        //                {
        //                    if (yyy > yh)
        //                    {
        //                        x = x1;
        //                        yyy = yh;
        //                    }
        //                    else
        //                    {
        //                        if (yyy < y0)
        //                        {
        //                            x0 = x;
        //                            yl = yyy;
        //                        }
        //                        else
        //                        {
        //                            x1 = x;
        //                            yh = yyy;
        //                        }
        //                    }
        //                }
        //                if (x == 1.0 | x == 0.0)
        //                {
        //                    mainlooppos = breaknewtcycle;
        //                    continue;
        //                }
        //                d = (aaa - 1.0) * Math.Log(x) + (bbb - 1.0) * Math.Log(1.0 - x) + lgm;
        //                if (d < Math.Log(_MIN_VALUE))
        //                {
        //                    break;
        //                }
        //                if (d > Math.Log(_MAX_VALUE))
        //                {
        //                    mainlooppos = breaknewtcycle;
        //                    continue;
        //                }
        //                d = Math.Exp(d);
        //                d = (yyy - y0) / d;
        //                xt = x - d;
        //                if (xt <= x0)
        //                {
        //                    yyy = (x - x0) / (x1 - x0);
        //                    xt = x0 + 0.5 * yyy * (x - x0);
        //                    if (xt <= 0.0)
        //                    {
        //                        mainlooppos = breaknewtcycle;
        //                        continue;
        //                    }
        //                }
        //                if (xt >= x1)
        //                {
        //                    yyy = (x1 - x) / (x1 - x0);
        //                    xt = x1 - 0.5 * yyy * (x1 - x);
        //                    if (xt >= 1.0)
        //                    {
        //                        mainlooppos = breaknewtcycle;
        //                        continue;
        //                    }
        //                }
        //                x = xt;
        //                if (Math.Abs(d / x) < 128.0 * _MACHINE_EPSILON)
        //                {
        //                    break;
        //                }
        //                i = i + 1;
        //                mainlooppos = newtcycle;
        //                continue;
        //            }
        //            else
        //            {
        //                mainlooppos = breaknewtcycle;
        //                continue;
        //            }
        //        }

        //        //
        //        // breaknewtcycle
        //        //
        //        if (mainlooppos == breaknewtcycle)
        //        {
        //            dithresh = 256.0 * _MACHINE_EPSILON;
        //            mainlooppos = ihalve;
        //            continue;
        //        }
        //    }

        //    //
        //    // done
        //    //
        //    if (rflg != 0)
        //    {
        //        if (x <= _MACHINE_EPSILON)
        //        {
        //            x = 1.0 - _MACHINE_EPSILON;
        //        }
        //        else
        //        {
        //            x = 1.0 - x;
        //        }
        //    }
        //    return x;
        //}

        /*************************************************************************
        Continued fraction expansion #1 for incomplete beta integral

        Cephes Math Library, Release 2.8:  June, 2000
        Copyright 1984, 1995, 2000 by Stephen L. Moshier
        *************************************************************************/
        private static double IncompleteBetaFractionExpansion(double a, double b, double x, double big, double biginv)
        {
            double xk = 0;
            double pk = 0;
            double pkm1 = 0;
            double pkm2 = 0;
            double qk = 0;
            double qkm1 = 0;
            double qkm2 = 0;
            double k1 = 0;
            double k2 = 0;
            double k3 = 0;
            double k4 = 0;
            double k5 = 0;
            double k6 = 0;
            double k7 = 0;
            double k8 = 0;
            double r = 0;
            double t = 0;
            double ans = 0;
            double thresh = 0;
            int n = 0;

            k1 = a;
            k2 = a + b;
            k3 = a;
            k4 = a + 1.0;
            k5 = 1.0;
            k6 = b - 1.0;
            k7 = k4;
            k8 = a + 2.0;
            pkm2 = 0.0;
            qkm2 = 1.0;
            pkm1 = 1.0;
            qkm1 = 1.0;
            ans = 1.0;
            r = 1.0;
            n = 0;
            thresh = 3.0 * MachineEpsilon;
            do
            {
                xk = -(x * k1 * k2 / (k3 * k4));
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;
                xk = x * k5 * k6 / (k7 * k8);
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;
                if (qk != 0)
                {
                    r = pk / qk;
                }
                if (r != 0)
                {
                    t = Math.Abs((ans - r) / r);
                    ans = r;
                }
                else
                {
                    t = 1.0;
                }
                if (t < thresh)
                {
                    break;
                }
                k1 = k1 + 1.0;
                k2 = k2 + 1.0;
                k3 = k3 + 2.0;
                k4 = k4 + 2.0;
                k5 = k5 + 1.0;
                k6 = k6 - 1.0;
                k7 = k7 + 2.0;
                k8 = k8 + 2.0;
                if (Math.Abs(qk) + Math.Abs(pk) > big)
                {
                    pkm2 = pkm2 * biginv;
                    pkm1 = pkm1 * biginv;
                    qkm2 = qkm2 * biginv;
                    qkm1 = qkm1 * biginv;
                }
                if (Math.Abs(qk) < biginv | Math.Abs(pk) < biginv)
                {
                    pkm2 = pkm2 * big;
                    pkm1 = pkm1 * big;
                    qkm2 = qkm2 * big;
                    qkm1 = qkm1 * big;
                }
                n = n + 1;
            }
            while (n != 300);
            return ans;
        }

        /*************************************************************************
        Continued fraction expansion #2
        for incomplete beta integral

        Cephes Math Library, Release 2.8:  June, 2000
        Copyright 1984, 1995, 2000 by Stephen L. Moshier
        *************************************************************************/
        private static double IncompleteBetaFractionExpansion2(double a, double b, double x, double big, double biginv)
        {
            double xk = 0;
            double pk = 0;
            double pkm1 = 0;
            double pkm2 = 0;
            double qk = 0;
            double qkm1 = 0;
            double qkm2 = 0;
            double k1 = 0;
            double k2 = 0;
            double k3 = 0;
            double k4 = 0;
            double k5 = 0;
            double k6 = 0;
            double k7 = 0;
            double k8 = 0;
            double r = 0;
            double t = 0;
            double ans = 0;
            double z = 0;
            double thresh = 0;
            int n = 0;

            k1 = a;
            k2 = b - 1.0;
            k3 = a;
            k4 = a + 1.0;
            k5 = 1.0;
            k6 = a + b;
            k7 = a + 1.0;
            k8 = a + 2.0;
            pkm2 = 0.0;
            qkm2 = 1.0;
            pkm1 = 1.0;
            qkm1 = 1.0;
            z = x / (1.0 - x);
            ans = 1.0;
            r = 1.0;
            n = 0;
            thresh = 3.0 * MachineEpsilon;
            do
            {
                xk = -(z * k1 * k2 / (k3 * k4));
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;
                xk = z * k5 * k6 / (k7 * k8);
                pk = pkm1 + pkm2 * xk;
                qk = qkm1 + qkm2 * xk;
                pkm2 = pkm1;
                pkm1 = pk;
                qkm2 = qkm1;
                qkm1 = qk;
                if (qk != 0)
                {
                    r = pk / qk;
                }
                if (r != 0)
                {
                    t = Math.Abs((ans - r) / r);
                    ans = r;
                }
                else
                {
                    t = 1.0;
                }
                if (t < thresh)
                {
                    break;
                }
                k1 = k1 + 1.0;
                k2 = k2 - 1.0;
                k3 = k3 + 2.0;
                k4 = k4 + 2.0;
                k5 = k5 + 1.0;
                k6 = k6 + 1.0;
                k7 = k7 + 2.0;
                k8 = k8 + 2.0;
                if (Math.Abs(qk) + Math.Abs(pk) > big)
                {
                    pkm2 = pkm2 * biginv;
                    pkm1 = pkm1 * biginv;
                    qkm2 = qkm2 * biginv;
                    qkm1 = qkm1 * biginv;
                }
                if (Math.Abs(qk) < biginv | Math.Abs(pk) < biginv)
                {
                    pkm2 = pkm2 * big;
                    pkm1 = pkm1 * big;
                    qkm2 = qkm2 * big;
                    qkm1 = qkm1 * big;
                }
                n = n + 1;
            }
            while (n != 300);
            return ans;
        }

        /*************************************************************************
        Power series for incomplete beta integral.
        Use when b*x is small and x not too close to 1.

        Cephes Math Library, Release 2.8:  June, 2000
        Copyright 1984, 1995, 2000 by Stephen L. Moshier
        *************************************************************************/
        private static double IncompleteBetaPowerSeries(double a, double b, double x, double maxgam)
        {
            double s = 0;
            double t = 0;
            double u = 0;
            double v = 0;
            double n = 0;
            double t1 = 0;
            double z = 0;
            double ai = 0;
            double sg = 0;

            ai = 1.0 / a;
            u = (1.0 - b) * x;
            v = u / (a + 1.0);
            t1 = v;
            t = u;
            n = 2.0;
            s = 0.0;
            z = MachineEpsilon * ai;
            while (Math.Abs(v) > z)
            {
                u = (n - b) * x / n;
                t = t * u;
                v = t / (a + n);
                s = s + v;
                n = n + 1.0;
            }
            s = s + t1;
            s = s + ai;
            u = a * Math.Log(x);
            if (a + b < maxgam & Math.Abs(u) < Math.Log(Double.MaxValue))
            {
                t = Gamma(a + b) / (Gamma(a) * Gamma(b));
                s = s * t * Math.Pow(x, a);
            }
            else
            {
                t = LogGamma(a + b, ref sg) - LogGamma(a, ref sg) - LogGamma(b, ref sg) + u + Math.Log(s);
                if (t < Math.Log(Double.MinValue))
                {
                    s = 0.0;
                }
                else
                {
                    s = Math.Exp(t);
                }
            }
            return s;
        }

        #endregion

        #region Error functions

        /// <summary>
        /// Error function
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double Erf(double x)
        {
            var s = Math.Sign(x);
            x = Math.Abs(x);

            if (x < 0.5)
            {
                var xsq = x * x;

                var p = 0.007547728033418631287834;
                p = 0.288805137207594084924010 + xsq * p;
                p = 14.3383842191748205576712 + xsq * p;
                p = 38.0140318123903008244444 + xsq * p;
                p = 3017.82788536507577809226 + xsq * p;
                p = 7404.07142710151470082064 + xsq * p;
                p = 80437.3630960840172832162 + xsq * p;

                var q = 0.0;
                q = 1.00000000000000000000000 + xsq * q;
                q = 38.0190713951939403753468 + xsq * q;
                q = 658.070155459240506326937 + xsq * q;
                q = 6379.60017324428279487120 + xsq * q;
                q = 34216.5257924628539769006 + xsq * q;
                q = 80437.3630960840172826266 + xsq * q;

                return s * 1.1283791670955125738961589031 * x * p / q;
            }

            if (x >= 10)
                return s;

            return s * (1 - ComplementErf(x));
        }

        /// <summary>
        /// Complementary error function
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public static double ComplementErf(double x)
        {
            if (x < 0)
                return 2 - ComplementErf(-x);

            if (x < 0.5)
                return 1.0 - ComplementErf(x);

            if (x >= 10)
                return 0;

            var p = 0.0;
            p = 0.5641877825507397413087057563 + x * p;
            p = 9.675807882987265400604202961 + x * p;
            p = 77.08161730368428609781633646 + x * p;
            p = 368.5196154710010637133875746 + x * p;
            p = 1143.262070703886173606073338 + x * p;
            p = 2320.439590251635247384768711 + x * p;
            p = 2898.0293292167655611275846 + x * p;
            p = 1826.3348842295112592168999 + x * p;
            
            var q = 1.0;
            q = 17.14980943627607849376131193 + x * q;
            q = 137.1255960500622202878443578 + x * q;
            q = 661.7361207107653469211984771 + x * q;
            q = 2094.384367789539593790281779 + x * q;
            q = 4429.612803883682726711528526 + x * q;
            q = 6089.5424232724435504633068 + x * q;
            q = 4958.82756472114071495438422 + x * q;
            q = 1826.3348842295112595576438 + x * q;

            return Math.Exp(-x * x) * p / q;
        }

        #endregion

        #region Bessel functions
        // from Accord.Net
        
        /// <summary>
        /// Returns the Bessel function of order 0 of the specified number.
        /// </summary>
        public static double BesselJ0(double x)
        {
            double ax;

            if ((ax = System.Math.Abs(x)) < 8.0)
            {
                double y = x * x;
                double ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
                    + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
                double ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
                    + y * (59272.64853 + y * (267.8532712 + y * 1.0))));

                return ans1 / ans2;

            }
            else
            {
                double z = 8.0 / ax;
                double y = z * z;
                double xx = ax - 0.785398164;
                double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
                double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                    + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                    - y * 0.934935152e-7)));

                return System.Math.Sqrt(0.636619772 / ax) *
                    (System.Math.Cos(xx) * ans1 - z * System.Math.Sin(xx) * ans2);
            }
        }

        /// <summary>
        /// Returns the Bessel function of order 1 of the specified number.
        /// </summary>
        public static double BesselJ(double x)
        {
            double ax;
            double y;
            double ans1, ans2;

            if ((ax = System.Math.Abs(x)) < 8.0)
            {
                y = x * x;
                ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1
                    + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
                ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74
                    + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
                return ans1 / ans2;
            }
            else
            {
                double z = 8.0 / ax;
                double xx = ax - 2.356194491;
                y = z * z;

                ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
                    + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
                ans2 = 0.04687499995 + y * (-0.2002690873e-3
                    + y * (0.8449199096e-5 + y * (-0.88228987e-6
                    + y * 0.105787412e-6)));
                double ans = System.Math.Sqrt(0.636619772 / ax) *
                    (System.Math.Cos(xx) * ans1 - z * System.Math.Sin(xx) * ans2);
                if (x < 0.0) ans = -ans;
                return ans;
            }
        }

        /// <summary>
        /// Returns the Bessel function of order n of the specified number.
        /// </summary>
        public static double BesselJ(int n, double x)
        {
            int j, m;
            double ax, bj, bjm, bjp, sum, tox, ans;
            bool jsum;

            double ACC = 40.0;
            double BIGNO = 1.0e+10;
            double BIGNI = 1.0e-10;

            if (n == 0) return BesselJ0(x);
            if (n == 1) return BesselJ(x);

            ax = System.Math.Abs(x);
            if (ax == 0.0) return 0.0;
            else if (ax > (double)n)
            {
                tox = 2.0 / ax;
                bjm = BesselJ0(ax);
                bj = BesselJ(ax);
                for (j = 1; j < n; j++)
                {
                    bjp = j * tox * bj - bjm;
                    bjm = bj;
                    bj = bjp;
                }
                ans = bj;
            }
            else
            {
                tox = 2.0 / ax;
                m = 2 * ((n + (int)System.Math.Sqrt(ACC * n)) / 2);
                jsum = false;
                bjp = ans = sum = 0.0;
                bj = 1.0;
                for (j = m; j > 0; j--)
                {
                    bjm = j * tox * bj - bjp;
                    bjp = bj;
                    bj = bjm;
                    if (System.Math.Abs(bj) > BIGNO)
                    {
                        bj *= BIGNI;
                        bjp *= BIGNI;
                        ans *= BIGNI;
                        sum *= BIGNI;
                    }
                    if (jsum) sum += bj;
                    jsum = !jsum;
                    if (j == n) ans = bjp;
                }
                sum = 2.0 * sum - bj;
                ans /= sum;
            }
            return x < 0.0 && n % 2 == 1 ? -ans : ans;
        }

        /// <summary>
        /// Returns the Bessel function of the second kind, of order 0 of the specified number.
        /// </summary>
        public static double BesselY0(double x)
        {
            if (x < 8.0)
            {
                double y = x * x;

                double ans1 = -2957821389.0 + y * (7062834065.0 + y * (-512359803.6
                    + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))));
                double ans2 = 40076544269.0 + y * (745249964.8 + y * (7189466.438
                    + y * (47447.26470 + y * (226.1030244 + y * 1.0))));

                return (ans1 / ans2) + 0.636619772 * BesselJ0(x) * System.Math.Log(x);
            }
            else
            {
                double z = 8.0 / x;
                double y = z * z;
                double xx = x - 0.785398164;

                double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
                double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                    + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                    + y * (-0.934945152e-7))));
                return System.Math.Sqrt(0.636619772 / x) *
                    (System.Math.Sin(xx) * ans1 + z * System.Math.Cos(xx) * ans2);
            }
        }

        /// <summary>
        /// Returns the Bessel function of the second kind, of order 1 of the specified number.
        /// </summary>
        public static double BesselY(double x)
        {
            if (x < 8.0)
            {
                double y = x * x;
                double ans1 = x * (-0.4900604943e13 + y * (0.1275274390e13
                    + y * (-0.5153438139e11 + y * (0.7349264551e9
                    + y * (-0.4237922726e7 + y * 0.8511937935e4)))));
                double ans2 = 0.2499580570e14 + y * (0.4244419664e12
                    + y * (0.3733650367e10 + y * (0.2245904002e8
                    + y * (0.1020426050e6 + y * (0.3549632885e3 + y)))));
                return (ans1 / ans2) + 0.636619772 * (BesselJ(x) * System.Math.Log(x) - 1.0 / x);
            }
            else
            {
                double z = 8.0 / x;
                double y = z * z;
                double xx = x - 2.356194491;
                double ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
                    + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
                double ans2 = 0.04687499995 + y * (-0.2002690873e-3
                    + y * (0.8449199096e-5 + y * (-0.88228987e-6
                    + y * 0.105787412e-6)));
                return System.Math.Sqrt(0.636619772 / x) *
                    (System.Math.Sin(xx) * ans1 + z * System.Math.Cos(xx) * ans2);
            }
        }

        /// <summary>
        /// Returns the Bessel function of the second kind, of order n of the specified number.
        /// </summary>
        public static double BesselY(int n, double x)
        {
            double by, bym, byp, tox;

            if (n == 0) return BesselY0(x);
            if (n == 1) return BesselY(x);

            tox = 2.0 / x;
            by = BesselY(x);
            bym = BesselY0(x);
            for (int j = 1; j < n; j++)
            {
                byp = j * tox * by - bym;
                bym = by;
                by = byp;
            }
            return by;
        }
        #endregion

        #region Temporary codes

        ///// <summary>
        ///// 
        ///// </summary>
        ///// <see cref="Sean Eddy's Squid"/>
        ///// <param name="x"></param>
        ///// <returns></returns>
        //public static double GammaLn(double x)
        //{
        //    if (x <= 0)
        //        throw new ArgumentOutOfRangeException(String.Format("Invalid x <= 0 in GammaLn(x = {0}).", x));

        //    int i;
        //    double xx, tx;
        //    double tmp, value;
        //    double[] cof = new double[11]{
        //        4.694580336184385e+04,
        //        -1.560605207784446e+05,
        //        2.065049568014106e+05,
        //        -1.388934775095388e+05,
        //        5.031796415085709e+04,
        //        -9.601592329182778e+03,
        //        8.785855930895250e+02,
        //        -3.155153906098611e+01,
        //        2.908143421162229e-01,
        //        -2.319827630494973e-04,
        //        1.251639670050933e-10
        //    };

        //    xx = x - 1.0;
        //    tx = tmp = xx + 11.0;
        //    value = 1.0;
        //    for (i = 10; i >= 0; i--) // sum least significant terms first
        //    {
        //        value += cof[i] / tmp;
        //        tmp   -= 1.0;
        //    }
        //    value = Math.Log(value);
        //    tx += 0.5;
        //    value += 0.918938533 + (xx + 0.5) * Math.Log(tx) - tx;

        //    return value;
        //}

        

        #endregion
    }
}
